% Producing Figure 8

close all
warning off

load('no_policy')

filename = 'rbc.mat';

temp      = load(filename,'solved_eta1_rhsee');
eta_rbc   = temp.solved_eta1_rhsee;
temp      = load(filename,'expmx');
expmx_rbc = temp.expmx;
temp      = load(filename,'npol');
npol_rbc  = temp.npol;
temp      = load(filename,'delta');
delta_rbc = temp.delta;

temp        = load(filename,'stss_y_rbc');
stss_y_rbc  = temp.stss_y_rbc;
temp        = load(filename,'stss_c_rbc');
stss_c_rbc  = temp.stss_c_rbc;
temp        = load(filename,'stss_h_rbc');
stss_h_rbc  = temp.stss_h_rbc;
temp        = load(filename,'stss_i_rbc');
stss_i_rbc  = temp.stss_i_rbc;
temp        = load(filename,'stss_k_rbc');
stss_k_rbc  = temp.stss_k_rbc;

%% Set up variables

options = optimset;

before = 4;
after  = 20;

stss_erk = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
             - (1-alpha)/(alpha*nu+1)*log(stss_k);
thre_z = ( stss_rkhats - stss_erk ) / sigma_rk;

shocks_run   = zeros(1,2+before+after);
shocks_norun = zeros(1,2+before+after);
shocks_run(2+before)   = thre_z-0.0001;
shocks_norun(2+before) = thre_z-0.0000;

%Exi = zeros(T+1);  % stochastic steady state

kv   = zeros(1,2+before+after);
rbv  = zeros(1,2+before+after);
nv   = zeros(1,2+before+after);
xiv  = zeros(1,2+before+after);
hv   = zeros(1,2+before+after);
cv   = zeros(1,2+before+after);
yv   = zeros(1,2+before+after);
iv   = zeros(1,2+before+after);
lv   = zeros(1,2+before+after);
pibv = zeros(1,2+before+after);
runv = zeros(1,2+before+after);
rhseev = zeros(1,2+before+after);
rkhat  = zeros(1,2+before+after);
rkhats = zeros(1,2+before+after);
prv  = zeros(1,2+before+after);
av   = zeros(1,2+before+after);
erk  = zeros(1,2+before+after);
z_temp = zeros(1,2+before+after);

x    = zeros(2+before+after,npol);
x_nonlog    = zeros(2+before+after,npol);
x_allnonlog = zeros(2+before+after,npol);

kv2   = zeros(1,2+before+after);
rbv2  = zeros(1,2+before+after);
nv2   = zeros(1,2+before+after);
xiv2  = zeros(1,2+before+after);
hv2   = zeros(1,2+before+after);
cv2   = zeros(1,2+before+after);
yv2   = zeros(1,2+before+after);
iv2   = zeros(1,2+before+after);
lv2   = zeros(1,2+before+after);
pibv2 = zeros(1,2+before+after);
runv2 = zeros(1,2+before+after);
rhseev2 = zeros(1,2+before+after);
rkhat2  = zeros(1,2+before+after);
rkhats2 = zeros(1,2+before+after);
prv2  = zeros(1,2+before+after);
av2   = zeros(1,2+before+after);
erk2  = zeros(1,2+before+after);
z_temp2 = zeros(1,2+before+after);

x2    = zeros(2+before+after,npol);
x_nonlog2    = zeros(2+before+after,npol);
x_allnonlog2 = zeros(2+before+after,npol);

kv3   = zeros(1,2+before+after);
hv3   = zeros(1,2+before+after);
cv3   = zeros(1,2+before+after);
yv3   = zeros(1,2+before+after);
iv3   = zeros(1,2+before+after);
rhseev3 = zeros(1,2+before+after);
rbv3    = zeros(1,2+before+after);
x3      = zeros(2+before+after,npol_rbc);

%% Simulation

kv(1)  = stss_k;
rbv(1) = stss_rb;
nv(1)  = stss_n;
lv(1)  = kv(1)/nv(1);
av(1)  = 0;

    for t=2:2+before+after
        av(t) = rho_a*av(t-1) + sigma_a*shocks_run(t);
        hv(t)  = ((1-alpha)/psi*exp(av(t))*kv(t-1)^(alpha))^(1/(alpha+1/nu));
        yv(t)  = exp(av(t)) * kv(t-1)^(alpha)*hv(t)^(1-alpha);

        rkhats(t) = log ( rbv(t-1) * (1-1/lv(t-1)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        rkhat(t)  = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                  + (1+nu)/(alpha*nu+1)*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t-1));
        xiv(t)  = exp(rkhat(t))/exp(rkhats(t));

        K_nonlog  = [log(kv(t-1)) rbv(t-1) log(nv(t-1)) xiv(t)];
        x0_nonlog = (K_nonlog).^expmx;
        x_nonlog(t,:)  = prod(x0_nonlog,2)';
        
        K_allnonlog  = [kv(t-1) rbv(t-1) nv(t-1) xiv(t)];
        x0_allnonlog = (K_allnonlog).^expmx;
        x_allnonlog(t,:)  = prod(x0_allnonlog,2)';

        K  = [kv(t-1) rbv(t-1) nv(t-1) xiv(t)];
        x0 = (log(K)).^expmx;
        x(t,:)  = prod(x0,2)';
        if xiv(t) >= 1
            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(lv(t-1)-1)*nv(t-1) - nv(t-1);

            if pibv(t) > 0
                rhseev(t) = exp(x(t,:)*solved_eta1_rhsee);
                rbv(t)    = exp(x(t,:)*solved_eta1_rbar);
                nv(t) = chi1*( chi0*pibv(t) + nv(t-1) ) + n0;
            else
                rhseev(t) = exp(x(t,:)*solved_eta3_rhsee);
                rbv(t)    = exp(x(t,:)*solved_eta3_rbar);
                nv(t) = chi1*(      pibv(t) + nv(t-1) ) + n0;
            end

        else
            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(1+lambda*(1-gamma))*(lv(t-1)-1)*nv(t-1) - nv(t-1);
            nv(t) = n_ss*(1 - nd/100);
            rhseev(t) = exp(x(t,:)*solved_eta2_rhsee);
            rbv(t)    = exp(x(t,:)*solved_eta2_rbar);
            runv(t) = 1;
        end

        cv(t) = rhseev(t)^(-1) + psi*hv(t)^(1+1/nu)/(1+1/nu);
        iv(t) = yv(t) - cv(t);
        kv(t) = (1-delta)*kv(t-1) + iv(t);

        lv(t)  = kv(t)/nv(t);

        erk(t) = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
               + (1+nu)/(alpha*nu+1)*rho_a*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t));

        rbv(t) = fzero('func_Rb',rbv(t),options,lv(t),erk(t),lambda,gamma,delta,sigma_rk);   % deposit interest rate

        rkhats_next = log ( rbv(t) * (1-1/lv(t)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        z_temp(t)     = (rkhats_next - erk(t)) / sigma_rk;
        cdf = normcdf(z_temp(t));
        prv(t) = cdf;

    end
    % End of T period simulation

%% Simulation, no run

kv2(1)  = kv(1);
rbv2(1) = rbv(1);
nv2(1)  = nv(1);
lv2(1)  = kv2(1)/nv2(1);
av2(1)  = 0;

    for t=2:2+before+after
        av2(t) = rho_a*av2(t-1) + sigma_a*shocks_norun(t);
        hv2(t)  = ((1-alpha)/psi*exp(av2(t))*kv2(t-1)^(alpha))^(1/(alpha+1/nu));
        yv2(t)  = exp(av2(t)) * kv2(t-1)^(alpha)*hv2(t)^(1-alpha);

        rkhats2(t) = log ( rbv2(t-1) * (1-1/lv2(t-1)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        rkhat2(t)  = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                  + (1+nu)/(alpha*nu+1)*av2(t) - (1-alpha)/(alpha*nu+1)*log(kv2(t-1));
        xiv2(t)  = exp(rkhat2(t))/exp(rkhats2(t));

        K_nonlog  = [log(kv2(t-1)) rbv2(t-1) log(nv2(t-1)) xiv2(t)];
        x0_nonlog = (K_nonlog).^expmx;
        x_nonlog2(t,:)  = prod(x0_nonlog,2)';
        
        K_allnonlog  = [kv2(t-1) rbv2(t-1) nv2(t-1) xiv2(t)];
        x0_allnonlog = (K_allnonlog).^expmx;
        x_allnonlog2(t,:)  = prod(x0_allnonlog,2)';

        K  = [kv2(t-1) rbv2(t-1) nv2(t-1) xiv2(t)];
        x0 = (log(K)).^expmx;
        x2(t,:)  = prod(x0,2)';
        if xiv2(t) >= 1
            pibv2(t)  = ( exp(rkhat2(t)) + 1-delta )*lv2(t-1)*nv2(t-1) - rbv2(t-1)*(lv2(t-1)-1)*nv2(t-1) - nv2(t-1);

            if pibv2(t) > 0
                rhseev2(t) = exp(x2(t,:)*solved_eta1_rhsee);
                rbv2(t)    = exp(x2(t,:)*solved_eta1_rbar);
                nv2(t) = chi1*( chi0*pibv2(t) + nv2(t-1) ) + n0;
            else
                rhseev2(t) = exp(x2(t,:)*solved_eta3_rhsee);
                rbv2(t)    = exp(x2(t,:)*solved_eta3_rbar);
                nv2(t) = chi1*(      pibv2(t) + nv2(t-1) ) + n0;
            end

        else
            pibv2(t)  = ( exp(rkhat2(t)) + 1-delta )*lv2(t-1)*nv2(t-1) - rbv2(t-1)*(1+lambda*(1-gamma))*(lv2(t-1)-1)*nv2(t-1) - nv2(t-1);
            nv2(t) = n_ss*(1 - nd/100);
            rhseev2(t) = exp(x2(t,:)*solved_eta2_rhsee);
            rbv2(t)    = exp(x2(t,:)*solved_eta2_rbar);
            runv2(t) = 1;
        end

        cv2(t) = rhseev2(t)^(-1) + psi*hv2(t)^(1+1/nu)/(1+1/nu);
        iv2(t) = yv2(t) - cv2(t);
        kv2(t) = (1-delta)*kv2(t-1) + iv2(t);

        lv2(t)  = kv2(t)/nv2(t);

        erk2(t) = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
               + (1+nu)/(alpha*nu+1)*rho_a*av2(t) - (1-alpha)/(alpha*nu+1)*log(kv2(t));

        rbv2(t) = fzero('func_Rb',rbv2(t),options,lv2(t),erk2(t),lambda,gamma,delta,sigma_rk);   % deposit interest rate

        rkhats_next = log ( rbv2(t) * (1-1/lv2(t)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        z_temp2(t)     = (rkhats_next - erk2(t)) / sigma_rk;
        cdf2 = normcdf(z_temp2(t));
        prv2(t) = cdf2;

    end

%% RBC model

kv3(1)  = stss_k_rbc;

    for t=2:2+before+after
        hv3(t)  = ((1-alpha)/psi*exp(av(t))*kv3(t-1)^(alpha))^(1/(alpha+1/nu));
        yv3(t)  = exp(av(t)) * kv3(t-1)^(alpha)*hv3(t)^(1-alpha);

        x0 = [log(kv3(t-1)) av(t)].^expmx_rbc;
        x3(t,:)  = prod(x0,2)';
        
        rhseev3(t) = exp(x3(t,:)*eta_rbc);

        cv3(t)  = rhseev3(t)^(-1) + psi*hv3(t)^(1+1/nu)/(1+1/nu);
        iv3(t)  = yv3(t) - cv3(t);
        kv3(t)  = (1-delta_rbc)*kv3(t-1) + iv3(t);
        rbv3(t) = 1 - delta + alpha*yv3(t)/kv3(t-1);

    end
    
%% plot resulsts

close all

figure('Name','consumption','Units','Inches','OuterPosition',[7 0 8 8])
plot(-before:after,cv(2:end)/stss_c*100-100,'Linewidth',4)
hold on
plot(-before:after,cv2(2:end)/stss_c*100-100,'r:','Linewidth',4)
hold on
plot(-before:after,cv3(2:end)/stss_c_rbc*100-100,'m-.','Linewidth',4)
hold on
plot(-before:after,zeros(1,1+before+after),'k:','Linewidth',2);
hold on
plot([0 0],[-20 5],'k-.','Linewidth',2);
set(gca,'FontSize',40);
xlim([-before after])
xticks(-before:4:after)
ylim([-20 5]);
yticks(-20:5:5)
grid on
legend({'run','no run','RBC'},'Fontsize',20,'NumColumns',3)
ytickformat('percentage')
%title('consumption','Fontsize',40)
xtickangle(0)

figure('Name','output','Units','Inches','OuterPosition',[7 0 8 8])
plot(-before:after,yv(2:end)/stss_y*100-100,'Linewidth',4)
hold on
plot(-before:after,yv2(2:end)/stss_y*100-100,'r:','Linewidth',4)
hold on
plot(-before:after,yv3(2:end)/stss_y_rbc*100-100,'m-.','Linewidth',4)
hold on
plot(-before:after,zeros(1,1+before+after),'k:','Linewidth',2);
hold on
plot([0 0],[-20 5],'k-.','Linewidth',2);
set(gca,'FontSize',40);
xlim([-before after])
xticks(-before:4:after)
ylim([-20 5]);
yticks(-20:5:5)
grid on
ytickformat('percentage')
%title('output','Fontsize',40)
xtickangle(0)

figure('Name','investment','Units','Inches','OuterPosition',[7 0 8 8])
plot(-before:after,iv(2:end)/stss_i*100-100,'Linewidth',4)
hold on
plot(-before:after,iv2(2:end)/stss_i*100-100,'r:','Linewidth',4)
hold on
plot(-before:after,iv3(2:end)/stss_i_rbc*100-100,'m-.','Linewidth',4)
hold on
plot(-before:after,zeros(1,1+before+after),'k:','Linewidth',2);
hold on
plot([0 0],[-20 5],'k-.','Linewidth',2);
set(gca,'FontSize',40);
xlim([-before after])
xticks(-before:4:after)
ylim([-20 5]);
yticks(-20:5:5)
grid on
ytickformat('percentage')
%title('output','Fontsize',40)
xtickangle(0)

figure('Name','labor','Units','Inches','OuterPosition',[7 0 8 8])
plot(-before:after,hv(2:end)/stss_h*100-100,'Linewidth',4)
hold on
plot(-before:after,hv2(2:end)/stss_h*100-100,'r:','Linewidth',4)
hold on
plot(-before:after,hv3(2:end)/stss_h_rbc*100-100,'m-.','Linewidth',4)
hold on
plot(-before:after,zeros(1,1+before+after),'k:','Linewidth',2);
hold on
plot([0 0],[-20 5],'k-.','Linewidth',2);
set(gca,'FontSize',40);
xlim([-before after])
xticks(-before:4:after)
ylim([-20 5]);
yticks(-20:5:5)
grid on
ytickformat('percentage')
%title('output','Fontsize',40)
xtickangle(0)
